Test divergence and curl module

This document present tests on divergence and curl module calculation using pygsf.

Preliminary settings

The modules to import for dealing with grids are:



In [1]:

    
from pygsf.mathematics.arrays import *



In [2]:

    
from pygsf.spatial.rasters.geotransform import *



In [3]:

    
from pygsf.spatial.rasters.fields import *

Divergence in 2D

The definition of divergence for our 2D case is:

\begin{align} divergence = \nabla \cdot \vec{\mathbf{v}} & = \frac{\partial{v_x}}{\partial x} + \frac{\partial{v_y}}{\partial y} \end{align}

Curl module in 2D

The definition of curl module in our 2D case is:

\begin{equation*} \nabla \times \vec{\mathbf{v}} = \begin{vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z} \\ {v_x} & {v_y} & 0 \end{vmatrix} \end{equation*}

so that the module of the curl is:

\begin{equation*} |curl| = \frac{\partial v_y}{\partial x} - \frac{\partial v_x}{\partial y} \end{equation*}

The implementation of the curl module calculation has been debugged using the code at [2] by Johnny Lin. Deviations from the expected theoretical values are the same for both implementations.

Vector field parameters: testing divergence

We calculate a theoretical, 2D vector field and check that the parameters calculated by pygsf is equal to the expected one.

We use a modified example from p. 67 in [3].

\begin{equation*} \vec{\mathbf{v}} = 0.0001 x y^3 \vec{\mathbf{i}} - 0.0002 x^2 y \vec{\mathbf{j}} + 0 \vec{\mathbf{k}} \end{equation*}

In order to create the two grids that represent the x- and the y-components, we therefore define the following two "transfer" functions from coordinates to z values:



In [4]:

    
def z_func_fx(x, y):

    return 0.0001 * x * y**3

def z_func_fy(x, y):

    return - 0.0002 * x**2 * y

The above functions define the value of the cells, using the given x and y geographic coordinates.

geotransform and grid definitions

Gridded field values are calculated for the theoretical source vector field x- and y- components using the provided number of rows and columns for the grid:



In [5]:

    
rows=100; cols=200



In [6]:

    
size_x = 10; size_y = 10



In [7]:

    
tlx = 500.0; tly = 250.0

Arrays components are defined in terms of indices i and j, so to transform array indices to geographical coordinates we use a geotransform. The one chosen is:



In [8]:

    
gt1 = GeoTransform(
    inTopLeftX=tlx, 
    inTopLeftY=tly, 
    inPixWidth=size_x, 
    inPixHeight=size_y)

Note that the chosen geotransform has no axis rotation, as is in the most part of cases with geographic grids.

vector field x-component



In [9]:

    
fx1 = array_from_function(
    row_num=rows, 
    col_num=cols, 
    geotransform=gt1, 
    z_transfer_func=z_func_fx)

vector field y-component



In [10]:

    
fy1 = array_from_function(
    row_num=rows, 
    col_num=cols, 
    geotransform=gt1, 
    z_transfer_func=z_func_fy)

theoretical divergence

the theoretical divergence transfer function is:



In [11]:

    
def z_func_div(x, y):
    
    return 0.0001 * y**3 - 0.0002 * x**2

The theoretical divergence field can be created using the function expressing the analytical derivatives z_func_div:



In [12]:

    
theor_div = array_from_function(
    row_num=rows, 
    col_num=cols, 
    geotransform=gt1, 
    z_transfer_func=z_func_div)

pygsf-estimated divergence

Divergence as resulting from pygsf calculation:



In [13]:

    
div = divergence(
    fld_x=fx1, 
    fld_y=fy1, 
    cell_size_x=size_x, 
    cell_size_y=size_y)

We check whether the theoretical and the estimated divergence fields are close:



In [14]:

    
assert np.allclose(theor_div, div)

Vector field parameters: testing curl module

We test another theoretical, 2D vector field, maintaining the same geotransform and other grid parameters as in the previous example. We use the field described in example 1 in [4]:

\begin{equation*} \vec{\mathbf{v}} = y \vec{\mathbf{i}} - x \vec{\mathbf{j}} + 0 \vec{\mathbf{k}} \end{equation*}

The "transfer" functions from coordinates to z values are:



In [15]:

    
def z_func_fx(x, y):

    return y

def z_func_fy(x, y):

    return - x

geotransform and grid definitions

Gridded field values are calculated for the theoretical source vector field x- and y- components using the provided number of rows and columns for the grid:



In [16]:

    
rows=200; cols=200



In [17]:

    
size_x = 10; size_y = 10



In [18]:

    
tlx = -1000.0; tly = 1000.0

Arrays components are defined in terms of indices i and j, so to transform array indices to geographical coordinates we use a geotransform. The one chosen is:



In [19]:

    
gt1 = GeoTransform(
    inTopLeftX=tlx, 
    inTopLeftY=tly, 
    inPixWidth=size_x, 
    inPixHeight=size_y)

Note that the chosen geotransform has no axis rotation, as is in the most part of cases with geographic grids.

vector field x-component



In [20]:

    
fx2 = array_from_function(
    row_num=rows, 
    col_num=cols, 
    geotransform=gt1, 
    z_transfer_func=z_func_fx)

vector field y-component



In [21]:

    
fy2 = array_from_function(
    row_num=rows, 
    col_num=cols, 
    geotransform=gt1, 
    z_transfer_func=z_func_fy)

theoretical curl module

The theoretical curl module is a constant value:

\begin{equation*} curl = -2 \end{equation*}

pygsf-estimated module of curl

The module of curl as resulting from pygsf calculation is:



In [22]:

    
curl_mod = curl_module(
    fld_x=fx2, 
    fld_y=fy2, 
    cell_size_x=size_x, 
    cell_size_y=size_y)

We check whether the theoretical and the estimated curl module fields are close:



In [23]:

    
assert np.allclose(-2.0, curl_mod)